﻿#region Using Statements
using System;
using System.Collections.Generic;
using System.Text;
using Microsoft.Xna.Framework;

#endregion

#region COPYRIGHT

/*
    Copyright (c) 2008, 2009, 2010 
        Roland Rosenkranz (Glatzemann@email.de)
*/

#endregion

#region LICENSE

/*
    This file is part of starLiGHT.Collision.

    starLiGHT.Collision is free software: you can redistribute it and/or modify
    it under the terms of the GNU Lesser General Public License as published by
    the Free Software Foundation, either version 3 of the License, or
    (at your option) any later version.

    starLiGHT.Collision is distributed in the hope that it will be useful,
    but WITHOUT ANY WARRANTY; without even the implied warranty of
    MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
    GNU Lesser General Public License for more details.

    You should have received a copy of the GNU Lesser General Public License
    along with starLiGHT.Collision.  If not, see <http://www.gnu.org/licenses/>.

    
    ADDITIONAL (commercial) LICENSES for starLiGHT.Collision are available on request.
*/

#endregion

namespace starLiGHT.Physic.TwoD
{
    public class SPHSimulation
    {
        internal float cellSpace;
        internal SmoothingKernel SKGeneral;
        internal SmoothingKernel SKPressure;
        internal SmoothingKernel SKViscos;
        internal float Viscosity;

        World world;
        List<FluidParticle> m_particles;

        public SPHSimulation(World world)
        {
            this.world = world;
            this.cellSpace = PhysicSettings.CELL_SPACE;
            this.Viscosity = PhysicSettings.VISCOSITY;
            SKGeneral = new SKPoly6(cellSpace);
            SKPressure = new SKSpiky(cellSpace);
            SKViscos = new SKViscosity(cellSpace);

            m_particles = new List<FluidParticle>();
        }

        public void addParticle(FluidParticle particle)
        {
            particle.ID = m_particles.Count;
            if (particle.neighbors == null) particle.neighbors = new Dictionary<int, FluidParticle>();
            particle.neighbors[particle.ID] = particle;
            m_particles.Add(particle);
        }

        public void removeParticle(FluidParticle particle)
        {
            m_particles.Remove(particle);
        }

        public int numParticles()
        {
            return m_particles.Count;
        }

        public List<FluidParticle> Particles
        {
            get { return m_particles; }
        }

        public void resetForces()
        {
            foreach (FluidParticle p in m_particles)
            {
                p.force = Vector2.Zero;
            }
        }

        public void update(Vector2 gravity, TimeStep step)
        {
            calculatePressureAndDensities();
            calculateForces();
            updateParticles(step, gravity);
            checkParticleDistance();
            foreach (FluidParticle p in m_particles)
            {
                Transform foo = new Transform(), bar = new Transform();
                p.Synchronize(world.m_contactManager.m_broadPhase, ref foo, ref bar);
            }
        }

        void calculatePressureAndDensities()
        {
            Vector2 dist;
            foreach (FluidParticle p in m_particles)
            {
                p.density = 0.0f;
                foreach (FluidParticle n in p.neighbors.Values)
                {
                    dist = p.position - n.position;
                    p.density += (p.mass * SKGeneral.calculate(dist));
                }
                p.updatePressure();
            }
        }

        void calculateForces()
        {
            Vector2 dist, force;
            float scalar;

            for (int i = 0; i < m_particles.Count; i++)
            {
                foreach (FluidParticle neighbor in m_particles[i].neighbors.Values)
                {
                    if (neighbor.ID < i)
                    {
                        if (neighbor.density > PhysicSettings.FLT_EPSILON)
                        {
                            dist = m_particles[i].position - neighbor.position;
                            // pressure
                            scalar = neighbor.mass * (m_particles[i].pressure + neighbor.pressure) / (2.0f * neighbor.density);

                            force = SKPressure.calculateGradient(dist);
                            force *= scalar;
                            m_particles[i].force -= force;
                            neighbor.force += force;

                            // viscosity
                            scalar = neighbor.mass * SKViscos.calculateLaplacian(dist) * Viscosity * 1f / neighbor.density;
                            force = neighbor.velocity - m_particles[i].velocity;
                            force *= scalar;
                            m_particles[i].force += force;
                            neighbor.force -= force;
                        }
                    }
                }
            }
        }

        void updateParticles(TimeStep step, Vector2 gravity)
        {
            foreach (FluidParticle particle in m_particles)
            {
                particle.update(step, gravity);
            }
        }

        void checkParticleDistance()
        {
            float minDist = 0.5f * cellSpace;
            float minDistSq = minDist * minDist;
            for (int i = 0; i < m_particles.Count; i++)
            {
                foreach (FluidParticle neighbor in m_particles[i].neighbors.Values)
                {
                    Vector2 dist = neighbor.position - m_particles[i].position;
                    float distLenSq = dist.LengthSquared();
                    if (distLenSq < minDistSq)
                    {
                        if (distLenSq > PhysicSettings.FLT_EPSILON)
                        {
                            float distLen = (float)System.Math.Sqrt(distLenSq);
                            dist *= (0.5f * (distLen - minDistSq) / distLen);
                            neighbor.position -= dist;
                            neighbor.positionOld -= dist;
                            m_particles[i].position += dist;
                            m_particles[i].positionOld += dist;
                        }
                        else
                        {
                            float diff = 0.5f * minDistSq;
                            neighbor.position.Y -= diff;
                            neighbor.positionOld.Y -= diff;
                            m_particles[i].position.Y += diff;
                            m_particles[i].positionOld.Y += diff;
                        }
                    }
                }
            }
        }
    }
}
